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Abstract 


Nomenclature 


A magnetic circuit model is developed for a PM motor 
for flywheel applications. A sample motor is designed 
and modeled. Motor configuration and selection of 
materials is discussed, and the choice of winding 
configuration is described. A magnetic circuit model is 
described, which includes the stator back iron, rotor 
yoke, permanent magnets, air gaps and the stator teeth. 
Iterative solution of this model yields flux linkages, 
back EMF, torque, power, and radial force at the rotor 
caused by eccentricity. Calculated radial forces are then 
used to determine motor negative stiffness. 


Introduction 


High speed flywheel systems are currently under 
development at NASA Glenn Research Center (GRC) in 
Cleveland, Ohio. Flywheels are being considered as an 
alternative to batteries and reaction wheels in space 
vehicles. Strengths of this technology include high 
energy density, long life, deep depth of discharge, and 
pulse power capability. Also, flywheels can be deployed 
in an array which provides both energy storage and 
attitude control. A system level flywheel testbed is 
operational at GRC. The system uses active magnetic 
bearings to provide a long life, low-loss suspension of 
the rotor. Permanent magnet (PM) motor-generators are 
used to store and retrieve energy from the rotor. 

This paper discusses a magnetic circuit model which 
will be used as a design tool for subsequent flywheel 
motor-generator designs. A sample motor-generator is 
developed to validate the design code. For simplicity, 
the motor-generator unit will be referred to simply as a 
motor for the remainder of the paper. 
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Design Approach 

A variety of methods are brought together to design the 
motor. The following is a quick summary of the design 
approach used: 

1 . Select motor configuration and materials 

2. Optimize fractional pitch 

3. Optimize winding distribution factor 

4. Define motor magnetic circuit model 

5. Solve model iteratively for flux (360°) 

6. Calculate flux linkages, and back EMF 

7. Calculate power and torque versus angle 

8. Calculate radial force versus rotor displacement 

Details describing implementation of each of the above 
steps are presented in the following sections. 


Motor Configuration 


densities (-2.2 Tesla), and for its high strength, which 
will allow high speed operation. 

The two candidate permanent magnet materials were 
samarium cobalt and neodymium iron boron; samarium 
cobalt was selected due to its superior high temperature 
performance. For the analysis performed in this paper, 
the magnets are assumed to have constant radial 
magnetization. 


Motor Winding Configuration 

Next, the fractional pitch is selected, to optimize 
harmonic content of the MMF harmonics. The pitch 
factor for a given harmonic is given as [1]: 



( 1 ) 


The motor configuration chosen for the flywheel design 
is a PM, radial gap, outer stator, inner rotor 
configuration; the motor will be a three phase, four pole 
design with 24 slots. This configuration was selected 
based on previous experience with motors used in 
flywheel applications; size of the motor is also defined 
by the specific application. The layout for this motor is 
shown schematically in Figure 1. 


All of the iron parts (the stator back iron, the stator teeth 
and the rotor yoke) are made of an iron cobalt alloy, 
which was selected for its ability to carry high flux 



The amplitude of the MMF harmonics is proportional to 

k pjj 

; we plot this versus the fractional pitch, which is 

h 

w 

defined as — , and select fractional pitch to optimize 
T P 

desired harmonic content. The plot is shown in 
Figure 2. Note that the harmonic number, h , has only 
odd values possible, since only odd harmonics exist in 
this motor. 

For this design, we wish to minimize third harmonic 
content. One full pitch for this motor is 24 slots/4 poles 
= 6 slots; thus reasonable options are 3/6, 4/6, 5/6 
fractional pitch or 6/6 full pitch. Since third harmonic 


Pitch Factor / h Vs. Factional pitch 



Figure 2. Pitch factor as a function of fractional pitch 
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content in Figure 2 is minimized with a pitch factor of 
0.67, a coil pitch of 4 was selected. 

The winding structure for the motor stator is double 
layer lap wound. The winding distribution factor is 
given as [ 1 ] : 


k dh 
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( 2 ) 


where zlz p is the fraction of a pole pitch that is occupied 
by the phase belt. With a four pole, three phase, 24 slot, 
4/6 fractional pitch machine, each 360 electrical degree 
phase belt has four coils, and the only design options 
available are a 60° or a 120° phase belt. Since 
considerably more voltage is generated by a 60° phase 
belt (k d i is 0.958 for a 60° phase belt versus 0.837 for a 
120° phase belt), the 60° phase belt was selected. 

With all three phases, 360 electrical degrees is spanned 
by 12 coils, at 30° increments. The electrical phase 
relationship of these coils is shown in Figure 3. 

Figure 4 shows the MMF’s of the coils when they are 
collected into 60° degree phase belts. Note that some of 
the coils are wound in with negative polarities. 

Finally, the winding layout of all three phases is shown 
in Figure 5. 



Figure 3. Coil electrical phase relationship 
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Figure 5. Winding layout for all three motor phases 


Magnetic Circuit Model 


Next, the motor is modeled as a magnetic circuit. 
Reluctance elements are included in the model for the 
back iron, stator teeth, air gap, and rotor yoke. MMF 
sources in the model include the PMs and the windings 
enclosing the stator teeth. The PM MMF is determined 
using the demagnetization curve for the magnetic 
material, and the windings that enclose each tooth are 
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Figure 6. Magnetic circuit model of motor 

summed together to determine the magnitude of the 
tooth MMF. The schematic of the motor magnetic 
circuit is shown in Figure 6. 


Using Eq. (5), define the following lumped reluctance 
quantities for the first 23 loops 

X n ~ —*^tooth(n) ~ ^ gap{n ) — 9? backiron(n ) 
~^tooth(n+ 1) ~^gap(n+ 1) — 9? yoke(n) 

Y n — ^tooth(n+ 1) ^gap(n+ 1) (^) 

= 91 tooth(n ) "I" 91 gap(n ) (9) 

Using Eq. (6), the lumped reluctance quantities for loop 
24 are 


This magnetic circuit is solved by mesh currents. 
Summing the MMF around the loop containing tooth 1 
and tooth 2 yields: 


X 2 4 ^tooth24 ^ gap24 9? backiron24 
~ 91 toothl ~ 91 gapl ~ ^ yoke 24 


( 10 ) 


^ 1 91 backiron 1 ^ tooth 2 F PM 2 

— (((q — 4>2 }(^tooth2 91 gap2 ) — *h9 ! yokel 
“Oh ~§24 toothl + ^gapl) + 

F PMl F toothl = 0 

Rearranging Eq. (3) yields: 

(~^backiron\ ~^tooth2 ~^gap2 
~ 91 yokel ~ 91 toothl ~ 91 gapl )*h 
+ tooth2 + ^gap2)§2 + (91 toothl + ^gap\)§24 
= F tooth2 + F PM 2 ~ F PM\ ~ F toothl 


X 2 4 91 tooth 1 9? gap 1 (11) 

Z 24 —%tooth24 ^gap24 (1^) 

Using Eq. (7) through (12), we now define the 
reluctance matrix as: 

"Xi Yj 0 0 Zj ' 

Z 2 X 2 Y 2 0 0 

0 Z 3 X 3 Y 3 0 ••• 0 (13) 

Y24 0 0 Z24 X24 


Using Eq. (4), a generalized expression for loops 2 to 
23 can be derived: 

( — backiron(n) ~^tooth(n+ 1 ) gap(n+ 1 ) ~^yoke(n) 

~^tooth(n) gap{n))§n 

i^tooth(n+ 1 ) + ^ gap(n+\)^§ n+\ (^) 

i^tooth(n) + ^ gap(n))§n-l 

~ F tooth(n+\) F PM{n+ 1) — F PM(n ) — F tooth(n) 


Now define the following MMF’s for the first 23 loops 
F n ~ F tooth(n+l) + F PM(n+ 1) ~ F PM(n) ~ F tooth{n) (14) 
and the loop 24 MMF is 

F 24 = ^W/*l + ^PM1 — F PM 24 ~ F tooth24 0 5 ) 

We now define the following two vectors 


Rewriting Eq. (5) for the case of loop 24 yields: 

( — ^backiron24 toothl ~ ^ gapl ~^yoke 24 

~ ^tooth2 4 ~ ^gap24 )§24 + (^toothl + 
(^tooth2 4 + 91 g ^24) ( t ) 23 
= F toothl + F PMl ~ F PM 24 ~ F tooth24 



(16) 
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F = 


F 1 

f 2 


F 24 


(17) 


Using Eq. (3) through (17), our model is described as 

9h|> = F (18) 

and by rearranging Eq. (18), solution for flux is thus 

<|) = <R _1 F (19) 

Because reluctances and MMF source values are 

dependant on flux density, this model is solved 
iteratively. First, an initial value of flux is assumed, and 
flux densities are calculated. Using the B-H curve for 
the metal parts (Figure 7) and knowledge of the 
mechanical gap size, reluctances for the metal parts and 
the air gap are calculated. The demagnetization curve 
(Figure 8) is used to generate the PM MMF at this flux 
density level, and the user-selected phase currents are 
used to calculate the tooth MMF. Equation 19 is then 
solved using the calculated source and reluctance 

values, providing new values of loop fluxes, and the 
process is repeated until the flux converges to a 
satisfactory limit. 


B-H Curve of Iron Cobalt Alloy 



Figure 7. B-H curve of selected iron-cobalt alloy 


Demagnitization Curve of Samarium Cobalt 



Figure 8. Demagnetization curve of PM material 


Simulation Results 


The motor simulation approach described in the 
previous section was used to calculate flux linkages and 
back EMF, motor torque and power, and rotor force. 
Details of these simulations are presented in the 
following sections. 

Flux Linkage and Back EMF 

The flux linkages are calculated by first running the 
simulation described in the previous section. Flux 
linking every individual coil is calculated, and 
individual coil fluxes in a phase belt are added together 
to obtain the flux linkages for each phase. A plot 
showing flux linkages for each phase is shown in 
Figure 9. 

By performing time derivatives of the flux linkages, 
back EMF for each phase is then calculated using 
Faraday’s law of induction: 


dk 

Vemf ~ dt 


( 20 ) 


The simulation results for the motor operating at the 
maximum flywheel speed of 60kRPM, displaying the 
back EMF for each phase, are shown in Figure 10. Note 
that, as desired, the third harmonic has been virtually 
eliminated. 
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Phase A, B, and C Back EMF’s at 60kRPM 



Motor Power and Torque 

The instantaneous power delivered to the back EMF can 
be expressed as follows: 

P ~ ia^a-emf z b v b-emf * c v c-emf (21) 

This power is equivalent to the mechanical power 
delivered to the rotor. 

P ~ Pmech ~ ®mech^ (22) 


where aw/* is the mechanical angular velocity of the 
rotor and x is the torque applied to the rotor. 
Combination of Eq. (21) and (22) yield the following 
expression for torque 

^ _ ia v a-emf h b v b-emf * c v c-emf (23) 

®mech 

In order to calculate torque as a function of electrical 
angle, the simulation was run under constant phase 
current, with phase A current set to 1.0 A (which 
automatically sets phase B and C currents at -0.5 A). 
Individual simulations were run over an entire electrical 
revolution at five degree increments, using an arbitrarily 
selected zero degree point. The resultant torque plot is 
shown in Figure 1 1 . Note that the maximum torque/amp 
generated by the motor is -0.018 Nm/amp, and that this 
maximum will be achieved when the motor controller is 
aligned at about 255 electrical degrees. The ripples in 
the torque plot are due to back EMF perturbations, 
which are caused by the stator teeth. 


Torque la=1, lb=— .5, lc=— .5 



Electrical Angle (Degrees) 

Figure 11. Generated torque versus electrical angle 

Next, we wish to calculate motor power. Given the 
selected slot size and the current density of copper, and 
assuming the slot is filled with copper at a packing 
factor of 0.5, the maximum phase current the motor can 
draw is 122 A. Based on this, the previous simulation 
was run once again with constant phase A current of 
122A (phase B and C at -61A). The result of this 
simulation, motor power as a function of electrical 
angle, is presented in Figure 12. Note that the maximum 
motor power under these conditions is -18KW. 
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Power la = 1 22, lb = -61 , Ic = -61 



Figure 12. Maximum power versus electrical angle 


electrical speed, which corresponds to harmonics that 
are at 24 times the mechanical spin speed (as expected 
with 24 teeth). Note also that the y component of force 
is zero; this means that all of the force on the rotor is in 
the radial direction. 

Next, the previous simulation was run seven more 
times, each with a ten mil rotor displacement from 

center, at 45, 90, 135 315 mechanical degree locations 

around the rotor radial center point. With each 
simulation, force vectors on the rotor were calculated 
for an entire revolution. Each of these points, along with 
the data from Figure 13, are plotted in Figure 14. Note 
that the ripple caused by the teeth is present in each of 
these data sets, which explains why each of the data sets 
appears as a line. Note also that the force on the rotor in 
each of these simulations is in the radial direction. 


Rotor Force Calculation 


The forces the stator produces on the rotor are equal 
and opposite to the forces on the stator tooth tips. The 
magnitude of the forces on the tooth tips is given by [2] : 


tooth-tip 


A tooth-tip 


d2 

n tooth-tip 


2-^0 


(24) 


Since there are 24 teeth in the rotor, the rotor 
mechanical angle for each of the tooth locations can be 
expressed as 


0 


n 


(n - 1) • 2 • 7i 
24 


(25) 


where n = 1.24. The general vector expression for force 
on the rotor at any of the tooth tips, from Eq. (24), is 

F n = 

A tooth-tip ' B tooth-tip r. . \ yi 

f ^ [i • sm(0 J+ j • cos(0„ )J 

2-po 

Force on the rotor caused only by the PMs is calculated 
next, for the case of a rotor which is slightly displaced 
from center. This simulation is run with zero phase 
currents, with the rotor displaced in the +x direction, for 
an entire electrical revolution, using Eq. (26). Note that 
although the air gap for the motor is 200 mils total, the 
actual rotor motion is limited to ±10 mils by the 
magnetic bearing system; thus, the displacement for this 
simulation was set to 10 mils. The result is presented in 
Figure 13. The x component of force has an average 
value of 13.45 lbs, with harmonics at 12 times the 



Electrical Angle (Degrees) 

Figure 13: Force versus angle with 10 mil 
x-direction displacement 


Force with Radial Dispacement of 10 Mils 
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Figure 14. x and y force components with 
radial offsets of 1 0 mils, at 8 locations 
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Next, force versus position plots were generated, for 
small and large displacement values. Once again, the 
simulation was run with no phase current; separate 
simulations were run as the rotor was moved in the +x 
direction. Figure 15 shows the force versus 
displacement plot for small displacements (less than 2 
mils) off center, and each plotted value is the average 
value of force over a 360 electrical degree rotation. This 
plot shows that the negative stiffness of the motor with 
the rotor close to the center is about 700 lbs/in. Note 
that for small displacements the negative stiffness is 
linear. Figure 16 shows force versus displacement for 
larger displacements. 

Force Vs. Displacement 

1 . 4 1 1 1 1 1 1 1 1 1 1 

1.2- - 


— 0 . 8 - 
.Q 

o 
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Figure 15. Force versus small displacement 


Force vs. Displacement 
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3 

O 

30 
20 - 
10 - 

1 1 1 1 1 
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Figure 16. Force versus large displacement 


Motor Design Summary 

A summary of the motor design described in this paper 
is provided in Table 1. 


Table 1. Motor design summary 
(dimensions in inches) 


Parameter 

Design value 

Poles 

4 

Turns/coil 

1 

Number of slots 

24 

Inner rotor diameter 

1.0 

Outer rotor diameter 

2.75 

Outer stator diameter 

6.50 

Magnet thickness 

0.15 

Mechanical air gap 

0.20 

Slot depth 

0.394 

Axial stack length 

1.50 

Yoke thickness 

0.725 

Back iron thickness 

1.281 

Slot pitch 

0.4123 

Tooth width at air gap 

0.2062 

Maximum speed 

60kRPM 

Motor power (kW 
@60kRPM) 

18 


Conclusions 

This paper describes a detailed motor-generator circuit 
model and simulation approach. A sample design and 
simulation are presented, and the simulation predicts 
motor back EMF, power, and torque versus electrical 
angle for the specified design. Additionally, force due 
to the permanent magnets (zero motor currents) is 
calculated versus rotor displacement, yielding negative 
stiffness relationships; for small displacements of the 
rotor (<2 mils), the negative stiffness was found to be 
linear. Also, simulations at 10 mil displacements all 
around the rotor center revealed that only radial forces 
were exerted on the rotor by the permanent magnets. 
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